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Nearly  all  flows  of  technological  importance  are  turbulent.  In  most  combustion 
devices,  turbulence  plays  a  vital  role  in  promoting  fluid  mixing  and  hence  accelerating 
combustion  rates.  While  experimental  techniques  are  making  significant  progress,  there  are 
unanswered  fundamental  questions  concerning  turbulent  mixing  that  remain  beyond  the 
reach  of  experimentalists.  These  unanswered  questions  are  an  obstacle  to  the  development 
and  testing  of  statistical  models  of  turbulence.  Answers  are  needed  so  that  these  models 
can  be  developed  into  reliable  design  tools. 

The  objective  of  this  work  is  to  answer  fundamental  questions  concerning  turbulent 
mixing,  through  numerical  simulations  to  turbulent  flows. 

The  details  of  the  mixing  process  are  difficult  to  study  experimentally  because  the 
smallest  length  and  time  scales  must  be  resolved.  Further,  it  is  not  sufficient  just  to 
measure  the  value  of  the  scalar  (e.g.  composition  or  temperature)  at  a  point:  rather  its  first 
two  derivatives  in  each  direction  are  the  prime  quantities  of  interest.  That  mixing  is  a 
microscale  process  compounds  experimental  difficulties,  but  it  provides  grounds  to  hope 
for  universality.  That  is  (according  to  classical  theory),  at  high  Reynolds  number,  the 
details  of  the  mixing  process  are  determined  by  a  few  macroscale  quantities,  but  are 
otherwise  independent  of  the  large  scales  of  the  flow. 

As  well  as  looking  at  a  scalar  and  its  derivatives  at  a  point,  mixing  and  combustion 
can  be  usefully  studied  in  terms  of  surfaces.  For  example,  in  a  simple  turbulent  diffusion 
flame,  reaction  is  concentrated  around  a  constant  property  surface.  In  premixed  turbulent 
combustion,  the  flame  can  form  a  highly  wrinkled  propagating  surface.  Surfaces  are 
difficult  to  study  experimentally  since  they  are  three-dimensional  and  encompass  all  scales. 


These  considerations  show  that  the  direct  numerical  simulation  of  turbulence 
provides  a  good  means  of  studying  turbulent  mixing.  All  scales  are  resolved,  and  any 
quantity  of  interest  can  be  determined  since  the  fields  are  known  for  all  positions  and  times. 
In  view  of  the  universality  of  the  small  scales,  it  is  sufficient  to  simulate  just  a  few  flows  - 
the  most  easily  computed  being  homogeneous  turbulence. 


As  is  described  in  the  following  sections,  we  have  performed  a  series  of  (64)3  and 
(128)3  direct  numerical  simulations  of  homogeneous  isotropic  turbulence.  From  the  data 
obtained  —  both  Eulerian  and  Lagrangian  —  we  have  been  able  to  address  several 
important  questions  in  mixing  and  combustion  theory. 

RESEARCH  OBJECTIVES 

1.  To  adapt  Rogallo's  direct  numerical  simulation  code  to  run  on  a  MASSCOMP 
minicomputer,  and  to  check  the  veracity  of  the  adaptions. 


2.  Similarly,  to  adapt  the  code  to  run  on  the  Cornell  National  Supercomputer  Facility 
machines. 


3. 


To  develop,  test  and  implement  a  numerical  algorithm  for  tracking  fluid  particles _ 

through  the  simulated  flow  field. 


4. 


To  develop,  test  and  implement  a  forcing  procedure  so  that  stationary  isotropic 
turbulence  can  be  studied.  :  r 
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5 .  To  study  the  mixing  of  a  passive  scalar,  primarily  through  an  examination  of  its  pdf. 

6.  To  study  the  Lagrangian  statistics  of  velocity,  acceleration  and  dissipation. 

7.  To  use  Lagrangian  time  series  of  the  velocity  gradient  to  address  mixing  and 
combustion  issues  in  terms  of  line,  surface  and  volume  elements. 

ACCOMPLISHMENTS 

1 .  We  have  translated  Rogallo's  code  into  FORTRAN  and  checked  the  results  in  detail 
against  Rogallo's  calculations.  The  MASSCOMP  has  a  built-in  array  processor 
capable  of  a  peak  rate  of  5  Mflops.  It  is  extremely  difficult  to  achieve  this  peak  rate, 
and  we  are  content  to  have  achieved  an  overall  rate  of  about  1  Mflops. 

At  one  stage  in  the  development  of  the  code,  the  i/o  was  taking  twice  as  long  as  the 
computation  time.  However,  with  the  aid  of  MASSCOMP  personnel,  we  incorporated 
synchronous  C  i/o  routines  that  have  virtually  eliminated  this  problem. 

On  this  dedicated  minicomputer  we  can  perform  (64)3  simulations  in  1  day. 

2.  The  Cornell  National  Supercomputer  Facility  operates  two  types  of  Supercomputer 
FPS264's  and  an  IBM  3090/600. 

In  the  second  year,  most  of  our  production  runs  were  performed  on  the  FPS264's. 
These  machines  have  sufficient  memory  for  us  to  perform  (16)3,  (32)3  and  (64)3 
calculations,  without  explicit  i/o.  To  perform  (128)3  calculations,  would  require  either 
explicit  i/o,  or  running  8  FPS264's  in  parallel.  Because  the  i/o  is  too  slow  we  did  not 
pursue  this  route. 

The  IBM3090/600  is  a  very  fast  scalar  machine  with  six  vector  processors  that  can  be 
run  in  parallel.  In  the  third  year  we  adapted  Rogallo's  code  to  this  machine,  by 
recoding  it  to  better  exploit  the  vector  and  parallel  features.  Using  six- way  parallelism, 
and  with  the  machine  being  lightly  loaded,  we  can  achieve  computing  speeds  in  excess 
of  those  on  the  Cray  XMP. 

With  this  latest  (and  hopefully  final)  version  of  the  code  we  have  performed  (128)3 
simulations. 

3 .  In  order  to  extract  Lagrangian  time  series  from  the  computed  flow  field,  a  fluid  particle 
tracking  algorithm  is  required.  One  has  been  developed  that  uses  a  second-order 
Runge-Kutta  scheme  in  time,  consistent  with  Rogallo's  code.  On  each  Runge-Kutta 
step,  interpolation  is  used  to  evaluate  the  fluid  particle  velocity  from  the  velocities  at 
the  grid  nodes  in  physical  space.  Initially  simple  (second-order  accurate)  linear 
interpolation  was  used.  But  numerical  tests  showed  that  this  is  insufficiently  accurate. 
We  have  developed  three  types  of  more  accurate  schemes:  Taylor  series  schemes; 
cubic  splines;  and,  schemes  in  which  the  interpolation  weights  are  optimized  with 
respect  to  the  turbulence  spectrum. 

A  comparative  testing  of  these  schemes  has  been  completed.  Such  testing  is  made 
quite  difficult  by  the  obvious  lack  of  an  analytic  solutions  for  turbulent  flow  that  can  be 
compared  directly  with  numerical  results.  We  have  performed  four  types  of  test: 
direct  measurement  of  the  interpolation  error,  determination  of  particle-tracking  errors 
in  a  simple  laminar  flow  for  which  an  analytical  solution  exists;  determination  of 


particle-tracking  errors  in  a  frozen  turbulence  field;  and,  determination  of  particle¬ 
tracking  errors  in  forced  isotropic  turbulence. 

The  results  show  that  the  cubic  spline  scheme  has  the  best  overall  performance.  Using 
it,  accurate  Lagrangian  time  series  can  be  extracted  without  unduly  increasing  the 
computational  cost  of  the  basic  Eulerian  code. 

A  paper  describing  the  particle-tracking  work  (Yeung  &  Pope  1987b)  has  been 
submitted  to  the  Journal  of  Computational  Physics. 

4.  Lagrangian  time  series  are  most  easily  analyzed  if  the  turbulence  is  stationary.  This 
can  be  achieved  by  forcing.  That  is,  energy  is  artificially  added  to  the  large-scale 
motions.  The  assumption  underlying  the  use  of  this  artifice  is  that  the  small-scale 
structures  in  the  turbulent  field  are  statistically  independent  of  the  large  scale  flow  and 
the  details  of  the  forcing,  being  confined  to  low  wavenumbers  and  large  scales,  do 
not,  therefore,  substantially  alter  the  values  of  the  small-scale  quantities.  This 
postulate  is  a  reasonable  one  for  high  Reynolds-number  flows  and  it  follows  from  the 
Kolmogorov  hypotheses.  But,  at  the  moderate  Reynolds-numbers  characteristic  of 
direct  numerical  simulations,  the  postulate  is  less  certain  since  the  energy-containing 
and  the  dissipation  regions  in  the  energy  spectrum  may  overlap. 

We  have  developed  a  forcing  scheme  which  is  different  from  the  previous  methods  of 
Siggia  and  of  Kerr.  Briefly  stated,  the  scheme  involves  adding  a  random  forcing  term 
to  die  momentum  equations.  The  Fourier  coefficients  of  this  term  are  Uhlenbeck- 
Omstein  processes  having  a  specified  time  scale  and  standard  deviation.  This  forcing 
is  restricted  to  the  low  wavenumber  modes.  The  energy  input  occurs  due  to  the 
correlation  of  the  Fourier  coefficients  of  the  forcing  and  velocity.  By  choosing  the 
time  scale  and  the  standard  deviation  of  the  forcing  we  can  control  the  time  scale  and 
average  magnitude  of  the  energy  input  » 

The  forcing  scheme  has  been  fully  tested  (Eswaran  &  Pope  1987a).  It  is  found  that 
statistically  stationary  turbulence  is,  indeed,  achieved.  A  method  of  choosing  the 
parameter  in  the  Uhlenbeck-Omstein  processes  has  been  established  that  allows 
turbulence  of  a  specified  Reynolds  number  to  be  generated.  Because  the  numerical 
algorithm  is  based  on  a  rectilinear  grid  in  physical  space,  some  anisotropy  is 
introduced.  Tests  show  that  this  is  small:  even  the  large-scale  motions  appear  closely 
isotropic.  Most  importantly,  it  is  found  that  the  details  of  the  forcing  do  not  affect  the 
small  scales.  Thus  conclusions  about  microscale  processes  drawn  from  studies  on 
forced  turbulence  have  general  applicability. 

The  paper  describing  this  study  (Eswaran  &  Pope  1987a)  has  been  accepted  for 
publication  in  Computers  and  Fluids. 

5.  Considering  the  importance  of  mixing  to  combustion,  it  is  remarkable  how  ignorant 
we  are  of  the  details  of  the  simplest  mixing  problem.  To  help  remedy  this  situation  we 
have  performed  Direct  Numerical  Simulations  of  a  passive  scalar  in  stationary 
turbulence. 

Consider  a  passive  scalar  field  <)>(&, t)  which  initially  is  composed  of  blobs  of  <p  =  -1 
and  an  equal  amount  of  blobs  of  <}>  =  1.  Initially,  then,  the  probability  density  function 

(pdf)  of  <j>  is  a  double-delta  function  distribution.  As  time  progresses,  due  to  mixing, 
the  delta  functions  become  smeared,  and  eventually  the  pdf  becomes  a  Gaussian 


4 


centered  on  <j>  =  0.  This  simple  mixing  problem  is  often  taken  as  a  test  of  mixing 
models  —  even  though  we  have  virtually  no  knowledge  of  how  the  pdf  evolves 
between  its  initial  and  final  stages. 

We  have  simulated  this  flow  and  extracted  joint  pdfs  of  $  and  its  dissipation  rate  e$, 
as  well  as  conditional  dissipation  rates.  The  results  are  described  in  a  paper  (Eswaran 
&  Pope  1987b)  which  has  been  accepted  for  publication  in  Physics  of  Fluids.  In 
accord  with  experimental  observations  we  find  that  the  decay  rate  of  the  scalar  depends 
strongly  on  the  initial  scalar  length  scale.  But  surprisingly,  it  is  found  that  the  shapes 
adopted  by  the  pdf  during  its  evolution  appear  to  be  independent  of  the  initial 
condition. 

These  data  will  prove  extremely  valuable  in  the  development  of  better  mixing  models. 

A  major  emphasis  of  our  work  is  to  extract  Lagrangian  statistics.  We  have  performed 
a  study  (Yeung  &  Pope  1987a)  in  which  Direct  Numerical  Simulations  are  used  to 
obtain  Lagrangian  time  series  of  velocity,  acceleration  and  dissipation.  From  these 
various  statistics  can  be  found. 

Our  numerical  simulations  have  been  successful  in  obtaining  a  wide-range  of  important 
Lagrangian  velocity  statistics,  many  of  which,  notably  the  cross-correlation  between 
dissipation  and  acceleration,  have  apparently  not  been  reported  before.  Numerical 
tests  show  that  on  643  grids,  the  resolution  is  adequate  to  obtain  accurate  second-order 
velocity  statistics  up  to  R\  at  least  SO.  However,  high-order  statistics  above  Rx  about 
30,  in  particular,  high-order  moments  of  velocity  increments,  may  be  underestimated. 
Nevertheless,  several  interesting  physical  conclusions  can  be  drawn. 

First,  a  Reynolds  number  of  Rx  =  SO  is  not  sufficient  high  for  the  Kolmogorov 
similarity  hypotheses  to  apply.  Second,  Lagrangian  velocity  increments  over  short 
time  intervals  are  observed  to  be  intermittent  and  highly  non  Gaussian,  the  more  so  as 
Reynolds  number  increase.  Third,  the  dissipation  is  found  to  be  close  to  having  a  log 
normal  distribution,  especially  for  higher  Reynolds  numbers.  Fourth,  dissipation  and 
acceleration  are  highly  correlated,  providing  evidence  of  the  close  connection  between 
dissipation  fluctuations  and  the  intermittency  exhibited  by  the  velocity  increments.  A 
transformation  is  time  normalized  by  the  local  Kolmogorov  scale  yields  more  nearly 
Gaussian  process. 

Several  important  topics  related  to  mixing  and  combustion  have  been  addressed  by 
studying  the  behavior  of  lines,  surfaces  and  volume  elements  in  turbulence.  This  is 
effected  by  extracting  time  series  of  the  Lagrangian  velocity  gradient  tensor. 

It  has  long  been  recognized  (at  least  since  Batchelor  1954)  that  mixing  is  closely 
related  to  the  deformation  of  infinitesimal  material  elements  —  lines,  surfaces  and 
volumes.  For  example,  the  stretching  of  material  surfaces  causes  neighboring  surfaces 
to  approach  each  other,  thus  increasing  gradients  and  mixing.  Indeed  in  some  cases 
(e.g.  early  times  or  large  Schmidt  numbers)  the  whole  mixing  problem  can  be  posed  in 
terms  of  the  behavior  of  material  surfaces. 

Through  DNS  we  have  determined  samples  of  the  deformation  tensor  which 
completely  describes  the  deformation  of  a  material  element  The  wealth  of  information 
obtained  is  described  by  Girimaji  and  Pope  (1988).  As  a  simple  example,  we  have 
studied  the  length  of  an  infinitesimal  line  element.  As  can  be  shown  analytically. 
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initially  the  length  does  not  increase,  because  there  is  no  alignment  between  the  line 
and  the  rate-of-strain  tensor.  But  eventually  the  length  increases  exponentially  with 
time. 

The  flamelet  model  of  turbulent  diffusion  flames  is  currently  quite  popular  even  though 
several  of  its  key  premises  remain  untested.  For  example,  it  is  assumed  that  scalar 

dissipation  %  is  uniquely  related  to  the  rate  of  strain  y  normal  to  the  stoichiometric 

surface;  and  further  it  is  assumed  that  there  is  but  a  small  probability  of  y  being 
positive  (corresponding  to  the  surface  being  compressed  in  its  own  plane).  We  have 
tested  these  hypotheses.  With  mild  assumptions,  a  stoichiometric  surface  is  well 
approximated  by  a  material  surface.  For  this  case  we  find  that  there  is  no  unique 

relationship  between  x  and  y,  and  that  y  is  positive  about  30%  of  the  time. 

A  useful  model  of  premixed  turbulent  flames  is  that,  at  any  instant,  the  flame  can  be 
regarded  as  a  surface  separating  reactants  from  products  (Pope  1987a,b,  Pope  & 
Cheng  1988).  This  surface  propagates  normal  to  itself  at  a  speed  w  relative  to  the 
fluid.  (To  a  first  approximation  w  is  the  laminar  flame  speed.)  When  w  is  small 
(compared  to  the  turbulence  intensity,  say)  it  has  been  assumed  that  the  propagating 
surface  remains  close  to  the  material  surface  on  which  it  originated.  This  is  certainly 
true  when  w*  (w  normalized  by  the  Kolmogorov  velocity  scale)  is  small;  and  then  z* 
(the  distance  between  the  surfaces,  normalized  by  the  Kolmogorov  length  scale)  is 
simply  proportional  to  w*.  From  DNS  we  have  computed  the  pdf  of  z*/w*.  From 
this  we  can  deduce,  for  example,  that  even  for  w*  =  1  there  is  a  significant  probability 
that  the  surfaces  will  separate  by  15  Kolmogorov  scales.  (Once  this  occurs,  the  linear 
analysis  used  breaks  down,  and  the  surfaces  are  likely  to  separate  yet  further.) 
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